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Abstract 

We measured the vertical pressure response function of a layer of sand submitted 
to a localized normal force at its surface. We found that this response profile depends 
on the way the layer has been prepared: all profiles show a single centered peak whose 
width scales with the thickness of the layer, but a dense packing gives a wider peak 
than a loose one. We calculate the prediction of isotropic elastic theory in presence 
of a bottom boundary and compare it to the data. We found that the theory gives 
the right scaling and the correct qualitative shape, but fails to really fit the data. 



PACS numbers : 
46.25.-y Static elasticity 
45.70.Cc Static sandpiles 
83.70.Fn Granular solids 



1 Introduction 

The statics of granular materials has been receiving recently a lot of attention, for a review 
see e.g. |T|]. An important issue is still to understand the mechanical status of an assembly 
of non-cohesive grains. In the small deformation limit, a classical viewpoint assumes a be- 
havior akin to an effective elastic medium. At a given confinement pressure, linear relations 
between stress and strain are measured and for larger strains, another picture is proposed 
based on a plastic modelling of the stress-strain relations. Therefore, for all practical pur- 
poses the available models used to describe granular matter in the quasi-static limit are of 
the elasto-plastic class with constitutive parameters determined empirically from standard 
triaxial tests ||. This elastic viewpoint is somehow corroborated by ultrasound propaga- 
tion experiments where, under large confining pressure, elastic moduli of p and s waves 
produced by a localized pulse can be measured ||. However, sound propagation measure- 
ments also evidence a strong 'speckle-like' component associated with the intricate contact 
force-paths topology or 'force chains' network. In a granular packing, contact forces of am- 
plitude larger than the average were found to organize in cells of sizes of about 10 grains 
diameters @] . The fragile character of these structures is even more obvious at low confin- 
ing pressure where, for example, ultra-small perturbations within the pile can completely 
modify the sound response spectrum ||. More generally, subtle self-organization prop- 
erties of the contact force network (also called the texture) were evidenced by thorough 
numerical studies by Radjai et al. J3J. 

At a macroscopic level, the pressure profile under the base of a sand heap built from a 
from a point source (i.e. from a hopper outlet), shows a minimum below the apex , but 
does not when the heap is constructed by successive horizontal layers || [7| . This surprising 
effect is currently viewed as a signature of the preparation history. Successive avalanches 
originated from the hopper outlet could have embedded a microscopic structure which is 
reflected macroscopically by an arching effect below the apex. 

The ability for a granular piling to change its texture (granular contact network, force 
chains geometry) in response to an external constraint, have cast legitimum suspicions on 
the fundamental validity of elasticity for packing of hard grains. For these reasons, a new 
class of models - called OSL for 'Oriented Stress Linearity' - was introduced by Bouchaud et 
al. H, which could explain remarkably well the sandpile data ||, as well as stress screening 
in silos ||10|| . These models have been the subject of a rather controversial debate ||11| . One 



of the reasons for that was the fact that they do not belong to the standard elasto-plastic 
class. As a matter of fact, they do not require the introduction of a displacement field, and 
the usual stress-strain relations are rather replaced by 'stress-only' ones which encode the 
history dependent state of equilibrium of the piling. In particular, the equations governing 
the stress distribution in these models are of hyperbolic type, which contrasts with the 
elliptic (or mixed elliptico-hyperbolic) equations of elastic (or elasto-plastic) modellings. 
An attractive feature of hyperbolic equations is that they have characteristic lines along 
which stress is transmitted, and which were argued to be the mathematical transcription 
of force chains that one can clearly see in granular systems JT^] . 



Measurements of the pressure response of a layer of sand submitted to a localized 



normal force at its surface soon appeared to be a way to discriminate between the different 
classes of models 0. Such a crucial experiment addresses at the deepest level, questions on 
the real mechanical status for a granular assembly. In elasticity, the shape of this pressure 
profile shows a single centered broad peak, whose width scales with the height h of the 
layer ||13|| . On the contrary, OSL models predict a response with two peaks (or a ring in 
three dimensions) on each sides of the overloaded point. Experiments [|14], |l|, [16] and 



simulations |L7| |18| have then been performed recently . Although the picture is far from 
being completely clear yet, the conclusions of these works can be roughly summarized 
as follows. For disordered systems, experiments definitively show elastic-like response, 
while regular packing exhibit OSL features. A third class of granular assemblies have also 
macroscopic equilibrium equations of the hyperbolic type: packings that can be prepared 
under the special isostaticity condition [|1^, [2(J, defining the uniqueness of the contact forces 



once the list of contacts is know |2^, |23fl . In pratice, this condition would correspond to a 
minimum number of contact per grains, such as frictionless contact forces for 2D random 
packing. A recent numerical result obtained for such an assembly explicitely shows an 
OSL-like propagation as a response to a localized force (at least on a scale up to 20 grains 
size) ||21|| . Note at last that, by contrast, the experiment presented in |24| rather claims a 



'diffusive' response function, in agreement with the stochastic scalar g-model |25| but these 
experiments were performed on a small size packing and in a rather specific geometry. 

In this paper, we present response function measurements obtained on large pilings 
made of natural sand. We show that it is possible to get rather different pressure profiles 
when preparing the packing with two different procedures. As in |14[], we found that elastic 



predictions give the right scaling and the correct qualitative shape but here, we perform 
a quantitative comparison between experimental data and isotropic elasticity predictions. 
We seek to answer precisely the question whether stress transmission properties can still 
be described using an isotropic elastic medium theory. What we found is that elasticity 
actually fails to really fit the data. 

The paper is organized as follows. In section [^ we expose how the measures are done 
and the way the data are calibrated. Then we show how different data-set can be obtained, 
depending on the sample preparation method. Section |3] is devoted to the calculation of the 
stress components at the bottom of an isotropic elastic layer of finite thickness h. In order 
to make this paper easier to read, the details of these calculations are given in appendices 
[A| and [B| for the two and three dimensional cases respectively. The comparison between 
the experiments and elasticity is done in section f|. At last, we conclude the paper with a 
discussion on the interpretation of the results in section |5|. 

2 The response function experiment 

The experimental technique that we use for the measure of the pressure response of a 
granular layer to the application of a localized vertical force F at its top surface has been 
described in details in ||14|| . The sketch of the set-up can be seen in figure [1| Briefly, the 



pressure P is measured by the tiny change of the electrical capacity of the probe due to the 
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Figure 1: Sketch of the experimental set-up. A localized vertical force F is applied on the top 
surface of the granular layer (z = 0). The corresponding pressure response on the bottom is 
measured at some distance r from that point. The vertical z-axis points downwards and we note 
h the thickness of the layer (between and 10 cm). We use natural 'Fontainebleau' sand whose 
typical diameter is ~ 0.3 mm. 



slight deformation of its top membrane. F is applied with a piston whose displacement is 
monitored and controlled to stay as small as possible (less than 500/im) . To gain sensitivity, 
F is modulated at a frequency /, and the probe signal is directed to a lock-in amplifier 
synchronized at / too. Any choice of / between 0.1 and 80 Hz gives the same result. 
We checked that the response P is linear in F. Both piston and probe have a surface in 
contact with the grains of ~ 1 cm 2 . The container is large enough (50 x 50 cm 2 ) to be 
able to neglect finite size effects due to the lateral walls. Its bottom plate is very rigid 
(Duraluminium, thickness 2 cm), and covered by a sheet of sand paper, in order to avoid 
sliding of the grains on the plate. 

We call r the horizontal distance between the piston and the probe. In order to measure 
the profile P(r), it is easier to vary r by moving the piston. In principle, the horizontal 
integral of the pressure profiles F* = J ( j~ OG dr2TirP(r) should be constant and equal to 
the force F applied at the surface. In fact, due to arching screening effects around the 
probe, this integral actually shows a large dispersion - see figure § - from an experiment 
to another and remains less than F. This screening effect is well known to be inherent 
to every mesurement of stresses in granular materials JTT|, but we could get rid of this 
problem of screening we use F* to renormalize the pressure measurements: P <— -^P. 
The accurate determination of F* is a crucial point when it comes to the quantitative 
comparison between experiments and theory. In particular, we were very carreful to take 
enough data points to have a good estimation of the experimental off-set at large r. 

A particular attention should be paid to the way the granular layer has been prepared. 
In order to observe different mechanical behavior, we chose two extreme procedures. There 
are schematized on figure |3|. The first one consists of making a packing as dense as possible. 
We add the sand by layers of 0.5 cm and after each layer, we push hard on the grains with 
a metallic plate. We can then reach a compacity of order of 0.7 - note that this 'layer 
by layer' procedure may create inhomogeneities in the density field. By contrast, to make 
the packing very loose, we first place a sieve on the bottom plate of the container, pour 
the grains into the box, and then gently pull up the sieve all through the grains. The 
corresponding compacity is of order of 0.6. 
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Figure 2: This plot shows the integral F* of the profiles P{r) taken from different experiments. 
The dispersion is quite large, but no systematic variation of F* with the layer thickness h can be 
evidenced. 



As already mentioned in |[L4|| , data taken from layers of several heights can be plotted 
together by rescaling lengths by h - which contradicts the 'diffusive' description as proposed 
by the g-model. The rescaled data h 2 P as a function of r/h can be seen on figure f|. 
This plot clearly shows that the response of the granular layer is 'history dependent': the 
pressure profile of a dense packing is much broader than that of a loose one. 

Bousinesq and Cerruti gave the expression of the stress response in the case of a isotropic 
semi-infinite elastic medium submitted to a localized and vertical unitary force F at r = 

. For the vertical pressure at point (r, z) , this expression is 



3F 



2tt (r 2 + z 2 )5/ 2 ' 



(1) 



This formula is independent of the Poisson coefficient v of the elastic material and thus does 
not have any adjustable parameter. It is therefore unable to reproduce the two different 
experimental pressure profiles. As a matter of fact, this function lays in between the two 
profiles - see figure (|. The pressure responses of a dense and a loose packing of sand are 
thus respectively broader and narrower than the standard elastic response profile. 

In the next section, we shall take into account the finite thickness of the layer and 
derive the corresponding expressions for the stresses, which will indeed depend on v. These 
expressions will then be quantitavely compared to the experimental data. 
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Figure 3: The packing of the grains has been prepared in two different ways: we either make it 
very dense (compacity ~ 0.7) by pushing hard on the grains with a metallic plate (a), or very 
loose (compacity ~ 0.6) by pulling up a sieve through them (b). 
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Figure 4: The pressure response profile depends on the way the system of grains was prepared: 
it is broader for a dense packing (empty circles) than for a loose one (filled circles). The response 
of a semi-infinite isotropic elastic medium (solid line) lays in between. For a given preparation, 
the experimental data collapse pretty well when renormalized by their F* factor - see text -, and 
rescaled by the height h of the layer. Here, we have plotted together measures on layers whose 
thickness varies from h ~ 30 to h ~ 60 mm. 



3 Elastic calculations 

In this section, we derive the expressions of the stress tensor components at the bottom of 
an isotropic elastic layer of finite thickness h. This calculation is not new, but as far as we 
know, the available litterature only provides numerical tables p7| that make fits difficult to 



perform. Such a calculation is a bit heavy, and we chose to present here its main lines only. 
The full details can be found in appendices |A| and [FJ for the two and three dimensional 
cases respectively. The formalism we use and the way the calculation is lead is directly 



inspired from |T3|, 28 



The stress state of an elastic material is described by its stress tensor components a^. 
At equilibrium, these quantities must verify the force balance equations: 

Vi(Tij = pgj (2) 

where p is the density of the material and gj the gravity vector. These relations are not 
enough to form a closed system of equations. An additional physical input is required. 
In plain elasticity theory, a displacement field w, is introduced - it measures the change 
in position, with respect to the reference state where no constrains are applied-, and 
the corresponding strain tensor Uij = \ (|^ + ^-J is related to the stresses via linear 
relations which involve two parameters which characterize this pure elastic material: its 
Young modulus Y and its Poisson coefficient v. 

It is possible to express all the equations in terms of the stress components only. Elim- 
inating the u^, one gets 

(lH-v)A^ + [l + (3-dH|^| = 0, (3) 

where d is the space dimension - these equations are not valid in the case of a non-uniform 
external body force. In particular, contracting i and j, we see that the trace of the stress 
tensor is an harmonic function, i.e. that Aakk = 0. These relations include (derivatives 
of) the force balance equations (Q). Taking the laplacian of (0), we also see that the a^ 
are bi-harmonic. 

The solutions of equations (^|) can be found in Fourier transforms. Let first focus on the 
two-dimensional (x, z) case. Because AAcr^ = 0, the general form of the vertical pressure 
o zz can be written as follows: 

r+oo r r -, r -, . 

a zz = j o dq cos(gx) { [A+(q) + qzB+(q)] e« z + [A~ z (q) + qzB^q)] e^} . (4) 

The expression for a xx is very similar. For the shear stress a xz , the cosinus factor should be 
replaced by sin(gx). In fact, only four of these twelve functions A's and 5's are independent. 
They are fully determined by the boundary conditions, and we can get this way explicit - 
but integral - expressions for the stresses. 

In elliptic problems like elasticity, stress or strain conditions must be specified on all 
the boundaries. Our aim here is to calculate the response of a layer of height h submitted 
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Figure 5: Stress response functions for a two-dimensional elastic material. The main plot com- 
pares the pressure profile of a semi-infinite system at depth z = h, with the response of a finite 
elastic layer of thickness h with either a rough or a smooth bottom. The first and third curves 
of each plot are independent of the Poisson coefficient v. For the second one, we chose v = 0.3 
but its shape depends only very weakly on the value of v - see inset where the maximum of the 
response has been plotted against v. The side plots show the other components of the stress 
tensor (shear and horizontal pressure). 



to a localized pressure at its top surface. We then suppose that the 'piston' which applies 
this overload is perfectly smooth and imposes, for example, a normalized (F = 1) gaussian 
profile Q(x) for the vertical pressure 



Q(x) 



V2na 2 



-x 2 /2a 2 



(5) 



where a is the adjustable width of this overload. The two conditions at the top are then 
(i) a zz (x, 0) = Q(x) and (ii) a xz (x, 0) = 0. Concerning the bottom, we assume that it is 
perfectly rigid, such that (iii) u z (x, h) = 0, and either very smooth or very rough. The last 
boundary condition is then (iv-a) a xz (x, h) = or (iv-b) u x (x, h) = respectively. 

Integrations in both smooth and rough bottom cases can be done numerically, and the 
corresponding results for the stress response a^ is plotted on figure || These integrations 
have been done for the specific choice of a very peaked gaussian overload: a = O.OOl/i - a 
quasi ^-function. For comparison, these plots are shown together with the Green's function 
of a vertically semi-infinite medium. 
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Figure 6: 3d equivalent of figure |5[ We chose again a = O.OOl/i and v = 0.3. The results are 
qualitatively the same as in two dimensions. 



All three a zz curves have rougly the same shape. This is not longer true when we plot 
the horizontal pressure instead of the vertical one: a xx is proportionnal to a zz (see equation 
(pl|)) for a rough bottom, but shows a double peak for a semi- infinite medium as well as for 
a smooth bottom - with a large negative central part. Negative values can be also seen for 
<y zz on figure ^| especially for the case of a smooth bottom. They are absolutely admissible 
for elastic material - no delamination between the material and the bottom is allowed. 

An interesting and rather non-intuitive point is that the finiteness of the elastic layer 
narrows the stress response. Only the response on the rough bottom depends on the value 
of the Poisson coefficient. We chose v = 0.3. This dependence is very weak for a zz - see 
inset of figure |5]. At last, it should be noted that all these curves scale with the height h. 

The axi-symmetric three-dimensional calculation is very similar, except that trigono- 
metric functions have to be replaced by Bessel ones in equations like (f|). Again, a numerical 
integration of the functions A's and i?'s can be done for a gaussian overload, and the cor- 
responding pressure profiles in both smooth and rough cases are plotted on figure |B|. They 
are also compared to the semi-infinite solution. The 3d results are qualitatively the same 
as in two dimensions. The wideness of the reponse function is non monotonic with the 
poisson ratio and presents a maximum for v ~ 0.27. A slight difference is that not only 
the 3d solution for the rough bottom depends on the Poisson coefficient, but the smooth 
bottom solution too. Again, this dependence is very weak for the vertical component of 
the stress tensor. 




4 A quantitative comparison 

In this section, we want to compare quantitatively the pressure response measurements with 
the elastic predictions. Among the two cases calculated in the previous section (rough and 
smooth bottom), the first one is the closest to our experimental situation - we checked that 
the shear stress at the bottom of the grain layer is finite. Therefore, only rough bottom 
elastic formulae are going to be used for the following fits. 

A set of experimental data is a file with three columns: the horizontal distance between 
the piston and the probe r^, the corresponding pressure measurement Pk and its typical 
dispersion AP&. There are N e ~ 15 such triplets for one pressure profile. As explained in 
section 0, the data have been renormalized by their factor F* in order to be of integral 
unity. 

We quantify the 'distance' between the experimental pressure profile P and the elastic 
predictions a zz by computing the average quadratic deviation E: 

(6) 

E = 1 would mean that a typical distance between theory and experimental data is one 
error bar. A value of E larger than 1 will then be considered as not good. Because a zz 
depends on the Poisson coefficient u, E is also a function of v. This function has a minimal 
value which gives the best fitting v. The precision of this value depends on the sharpness 
of this minimum. 

The results of our fits are gathered together in figure |9|. The results can be summarized 
as follows. In the case of a dense packing, the experimental response function is too wide to 
be well fitted by an elastic curve - see figure [7|. As a matter of fact, we get typically E ~ 4 
for the best fit. The corresponding Poisson coefficient value does not then have any real 
meaning. For the loose packing, the situation is different. The experimental data, though 
closer to the elastic response, lay on a curve which is too narrow to be properly fitted. The 
best value of v is then v = | which gives the most narrow elastic response, corresponding 
to E ~ 2. Interestingly however, if one allows v to exceed the standard limit v — ~, one can 
fit the data pretty well - see figure || This exceeding can be done mathematically because 
the qualitative shape of the stress profiles calculated with the elasticity theory changes for 
v > | only, leading beyond this value to oscillatory behaviours - see appendices. Note that 
the dilatancy effect in granular material has sometimes been argued to be somehow encoded 
by a Poisson coefficient larger than |, which is the 'incompressibility' limit. Although we 
do not think that dilatance can be treated with the concepts of reversible elasticity, our 
results would contradict such an argument because only dense packings dilate, loose ones 
on the contrary contract. 

As we said in section |2|, the piston which applies the overload at the top surface of the 
layer, as well as the pressure probe have an area of ~ 1 cm 2 . Taking into account the finite 
size of the piston for the fits was easy in our elastic formalism. Indeed, we assumed that 
the overload had a gaussian profile of adjustable width a, which was then simply set to 
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Figure 7: Fit of the data obtained on a dense packed granular layer. It is rather poor because 
the elastic response cannot get wide enough, the inset shows the deviation E versus v. v = 0.27 
corresponds to E = 3.8. 
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Figure 8: Fit of the loose packing data. The best Poisson coefficient exceeds the usual v 
limit, v = 0.58 corresponds to E = 0.6, but E = 1.5 when v = 0.50. 
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Figure 9: This figure shows (a) the best elasticity fit quality E and (b) the corresponding value 
for the best poisson ration v as a function of the compacity of the packing. Above the dash line 
E = 1, the fits are not in good agreement with the data. Only point under the dash line v = \ 
are in principle valid in the elastic framework. 



the piston diameter. By contrast, taking into account that of the pressure probe requires 
the convolution of the elastic formulae with a disk of finite diameter. This calculation is 
much less easy and makes the computation of the stress profiles difficult to do. In fact, we 
checked on few data sets that these two finite size effects are not very important as soon 
as the layer thickness h is larger than ~ 30 mm, i.e. 3 piston/probe diameters. The actual 
values of E and v that come out from these modified fits are a bit different - slightly better 
- than that of figure || but the conclusions written in the previous paragraph keep exactly 
the same. 



5 Discussion and conclusions 

As it was shown that the stress profile under a sandpile does or does not has a 'dip' below 
the apex of the pile, we found that the pressure response of a layer of sand submitted to 
a localized normal force at its top surface depends on the way this layer has been built - 
its 'history'. The response is rather wide when the grain packing is made very dense and 
compact, but it is more narrow when the layer is loose. For a given height, the maximal 
value of the pressure is approximately twice smaller in the first than in the second case. 
The predictions of isotropic elasticity theory, even when taking into account the finite 
layer thickness, agree poorly with the experimental data. However note the puzzling result 
that, taking a Poisson coefficient larger than the usual limit |, can fit rather well the 
experimental data for the loose packing preparation. We have no interpretation of this 
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fact besides concluding for the non adequacy of the isotropic elastic picture for piling 
prepared using 'dense' or 'loose' filling procedures. 

Although we present all our results in terms of dense or loose packing, the compacity in 
itself is certainly not a good control parameter. Rather, a natural way to improve the fits 
within this elastic framework would be to take into account a possible anisotropy of the 
material. A classical example is the so-called 'aelotropy' which occurs when the vertical 
symmetry axis has different mechanical properties than the horizontal directions. Such 
an anisotropy has five independent parameters: two Poisson coefficients (a vertical and a 
horizontal one), two Young moduli (idem) and a shear modulus. A theory with so many 
parameters will for sure fit our data. Indeed, the shear modulus has been found to be of 
strong influence on the shape of the response function p9f . 

We think however that a standard elastic description of granular materials is unsat- 
isfactory: a proper definition of the kinetic variables may be problematic for systems of 
hard particles. As a matter of fact, the link between the local microscopic movement of 
the grains and the possible corresponding large scale displacement field is a current sub- 
ject of research [J3UJ. Recent numerical simulations of frictionless disks even suggest that 



the stress-strain relation might not converge to a well defined curve for larger and larger 
systems [ J3lH . 

A rather striking feature of granular systems is the presence of force chains. These 
chains support most of the weight of the grains, and their geometrical characteristics - 
length, orientation - is the signature of the history of the system. In a recent paper [K2| , 
some of us with others have shown that it is possible to get pseudo-elastic equations from 
a simple model of - perfectly rigid - force chains which can split or merge at some 'defects' 
of the grain packing. In this model however, the stress tensor as well as the vector field 
which plays the role of the displacement in elasticity can be both built from the angular 
distribution of force chains. The specification of the boundary conditions is therefore a 
non trivial issue on which we are currently working. There are two main advantages 
in this new approach. First, no real displacement field is needed, and second, it allows 
to calculate the pseudo-elastic coefficient from microscopic quantities - the force chains 
angular distribution. The idea is then to introduce some anisotropy in this distribution, 
and see which kind of anisotropic pseudo-elastic equations we get out of it - note that this 
work would be very close in spirit to e.g. [[33], |34| were they try to link local geometrical 



variables (orientation of contacts between grains) with the mechanical properties at larger 
scales. The fit of our data would then give an information on the local structure of the 
packing. More response function experiments are thus planed to be performed with new 
preparation history, in particular with shearing or avalanching procedures in order to be 
able to come back to the yet unresolved sandpile 'dip' problem. 
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A Appendix: the two-dimensional elastic calculation 

For a two-dimensional elastic layer, we have three independent stress tensor components: 
the pressures a xx and a zz , and the shear a xz . x is the horizontal axis, and z is the vertical 
one, pointed downwards. The continuity equations (Q) can then be explicitly written down 
as follows: 



d z a zz + d x a xz = pg 
d z cr xz + d x a xx = 0. 



(7) 



Besides these two equilibrium equations, an additionnal and independent equation is re- 
quired to solve the problem. Among equations (§), the simplest is 



A(a zz + a xx ) = 0. 
It is natural in this context to introduce the new variables 



(9) 



T 

T 

D 



O ZZ T & xx 
Oxz 

a,-. — a n 



' ZZ 



'xx- 



Using equations (|7|) and (|), it is easy to show that these new functions verify 



AT 

At 
d x d z D 





-d x d z T 

(d 2 z ~ d 2 x )r. 



(10) 

(11) 
(12) 



(13) 
(14) 
(15) 



A standard mathematical base of harmonic functions is the product of trigonometric 
functions with exponentials. We shall keep to x <-> — x symmetrical situations such that 
we can look for a solution of the type 



T 

T 

D 



B x z + Ci+ dq cos(qx) a(q)e qz + b(q)e qz 
Jo L 

1 r+oo 



B 2 x + dq sm(qx) < c(q)e gz + d(q)e qz 
2 Jo I L 

2pgz - (B 1 + B 2 )z + C 2 



+ 2 qZ 



a(q)e qz + b(q)e 



qz 



(16) 
(17) 



- / dq cos{qx) {qz [a{q)e qz - b(q)e~ qz ] + 2 [c{q)e qz - d{q)e' qz ] } (18) 

where the constants Bi, B 2 , C\ and C 2 , as well as the the functions a(q), b(q), c(q) and 
d(q) are to be determined by the boundary conditions. 
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Boundary conditions 

We suppose that the top surface is submitted to a localized vertical and unitary overload, 
for example a normalized gaussian profile Q(x) for the vertical pressure: 



Q(x) 



1 



-x 2 /2a 2 



V2 



ITU* 



(19) 



o is the adjustable width of this overload. The two conditions at the top are then (i) 
&zz{x, 0) = Q(x) and (ii) a xz (x,0) = 0. Concerning the bottom, we assume that it is 
perfectly rigid, such that (iii) u z (x, h) = 0, and either very smooth or very rough. The last 
boundary condition is then (iv-a) a xz (x, h) = or (iv-b) u x (x, h) = respectively. 

It is convenient for our calculation to transform the displacement conditions (iii) and 
(iv-b) into stress conditions. For that purpose, one can take derivatives of condition (iii) 
with respect to x and introduce stress components via the strain-stress relations. Using 
also equilibrium equations (0) and (||), one finally gets for the condition (iii): 



(2 + v)d x a xz {x, h) = d z a xx (x, h) - vpg. 
A similar calculation leads to the following new condition (iv-b): 

(x,h) = ua zz (x,h). 

and 



o. 



Since we look for a solution of the form of (|16D , flD,, 
the function s(q) such that 

r+oo 

Q(x) = I dq cos(qx)s(q), 



which, for the specific gaussian choice ( |19"D leads to 

1 



s(q) 



-e- a2q2/2 . 

71 



(20) 

(21) 
it is natural to introduce 

(22) 
(23) 



Solution for a smooth bottom 

We switch gravity off since it is not of interest for the calculation of the response function. 
A simple term to term identification in boundary conditions (i), (ii), (iii) and (iv-a) leads to 
vanishing coefficients Bi and Cj, and to the four following linear equations for the unknown 
functions a(q), b(q), c(q) and d(q): 



~[a(q)+b(q)]-c(q) + d(q) 
c(q) + d(q) 



:i+v) 



c(q) + -qha(q) 
c(q) + -qha(q) 



e qh + 



e qh + 



d(q) + -qhb(q) 
d{q) + -qhb(q) 



-qh 



- qh 



s{q) (24) 

(25) 

a{q)e qh - b{q)e- qh (26) 



0. 



(27) 
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whose solution is 



o(g) = 2s(q)- 



sinh(qh)e 



-qh 



smh(2qh) + 2qh 
smh(qh)e qh 



sinh(2qh) + 2qh 
c(q) = -d(q) = --s(q) — 



2 v smh(2qh) + 2qh' 



(28) 
(29) 
(30) 



Putting these relations back to the equations (|T6|) , (|T7j) and fllBD , we get explicit - integral - 
expressions for the stress components. Note that the function a, b, c and d are independent 
of v. 

Solution for a rough bottom 

For the case of a rough bottom, the condition (iv-a) has to be replaced by the condition 
(iv-b). It means that the last equation of the system (|24| - |2"7D has to be changed into 



c(q) + -qha(q) 



e qh + 



d(q) + -qhb(q) 



-qh 



11-1/ 



a(q)e qh + b(q)e- qh ] (31) 



2 1 + 1/ 

and the resolution of these four linear equations leads this time to the following solution 

f-(q)+2qh 



a(q) 
b(q) 
c(q) 



2s(q) 
2s(q) 
-d(q) 



f + (q)f4q)+M 2 h 2 

U(q)-2gh 
U( q )Uq)+^ 2 h 2 



s{q) 



f+(q) + f-(q) 

U(q)f-(q) + ±q 2 h 2 



where the functions /+ and /_ are defined by 

/±(g) = i + 



3-i/ 



±2qh 



(32) 
(33) 
(34) 

(35) 



This time, there is a dependance in v. In principle, the Poisson coefficient should be less 
that unity in 2d. However, these functions really change behaviour only for v > 3, leading 
to oscillatory stresses - they however develop negative parts as v — > 3. 



Semi-infinite medium 

For comparison, in the case of a semi-infinite medium submitted to a ponctual and vertical 
unitary force at x = 0, the stress components are given |13| by 



O", 



7T (X 2 + Z 2 ) 2 
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(36) 



zx 



o. 



<?<, 



7T (x 2 + Z 2 ) 



2\2 



XZ 



7T (x 2 + Z 2 ) 



T\2' 



(37) 
(38) 



B Appendix: the three-dimensional case 



The calculation in the three-dimensional case is very similar to the 2d one. We shall keep to 
axi-symmetric situations so that the stress tensor has only four non-zero components: the 
pressures a zz , <r rr , a 00 and the shear a rz . z is again the vertical axis pointing downwards, 
and (r, 9) are the horizontal planar coordinates. 

The equations we want to solve are simpler with the new functions 

(39) 
(40) 
(41) 
(42) 

which must verify 

(43) 
(44) 

(45) 

(46) 

The two last equations are the explicit forms of the force balance equations (0), and we 
got the two first ones from equations (|3|). 

The corresponding general solutions involve Bessel functions of the first kind Jo, J\ and 
J 2 , and read: 



T = 


0" 2 z + C"Vr + O"06» 


T = 


O rz 


S = 


a rr + a 00 


D = 


a rr — a 00, 


AT 


= 


(l + v)AS 


= dlT 


d r r + - 
r 


= -d z (T-S)+pg 


d r D + 2— 


= -d r S-2d z T. 



T 

S 



r+00 r 

T x z + T 2 + dq J {qr) \a{q)e qz + b{q)e' qz 
Jo L 

Siz + S 2 

+ J% J {qr) J [c{q)e qz + d{q) e - qz ' ' ' 



+ 



2 1 + v 



qz 



a(q)e qz - b(q)e 



-qz 



(47) 
(48) 



+00 



I 1 u(z) 
-pgr + -r(S 1 -T 1 ) + ^+ I dqJ 1 (qr)[[c(q)-a(q)]e qz + [d(q)-b(q)}e- 

II r Jo L 



qz 



1 



D 



2 1 + isJo 
v(z) du(z) 



dq Ji(qr) a(q) [1 + qz]e qz — b(q) [1 — qz]e~ 



qz 



d + J q dq J 2 (qr) [[2a(q) - c(q)]e qz + [26(g) - d(q)] e - qz 



1 



1 + v Jo 



dqJ 2 (qr) 



a{q)[2 + -\qz]e qz + b{q)[2- l -qz}e- qz 



(49) 



(50) 
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where the constants 7\, T 2 , Si and S2, as well as the functions u(z), v(z), a(q), b(q), c(q) 
and d(q) are, again, to be determined by the boundary conditions. 

Boundary conditions 

As in the 2d case, we want to impose at the surface (i) an overload cr zz (r, 0) = Q{r), but (ii) 
no shear a rz (r, 0) = 0. At the bottom, the vertical displacement must vanish u z (r, h) = 
(iii), and we shall study the two cases, very smooth a rz (r,h) = (iv-a) or very rough 
u r (r, h) = (iv-b) bottom. 

The 3d equivalent of (|19D is now 



Q(r) = 7^- 2 e- r2/2 ° 2 . (51) 

Looking at the general form of the solution (|47| - |50D , it is natural to introduce the function 

s(q) defined by the relation 

r+00 

Q{r)= / dqJ (qr)s(q) (52) 

Jo 

which, for Q given by fl5"T|) gives 

s(q) = ^qe~^ 2/2 . (53) 

Taking derivatives of the conditions on displacements and using stress-strain relations, 
one can again transform these conditions into relations between stress components only. 
It is easy to show that the condition (iii) can be written as 

2(1 + u)d r r = - (1 + v)d z {S + D)- ud z T, (54) 

and that the condition (iv-b) gives 

(1 + u)(S + D) = 2uT (55) 

(these last two relations are only valid at z = h). 

Solution for a smooth bottom 

Again, as we are interested to the response of the elastic layer to this overload, gravity 
is switched off. The four conditions (i)-(iv-a) then give four equations for the unknown 
functions a(q), b(q), c(q) and d(q), all other functions and constants being zero. These 
equations are: 

a(q) + b(q) - c(q) - d(q) = s(q) (56) 

c(q)-d(q) = ^±^ [ fl (g) - 6(g)] (57) 

2(1 + v) [c{q)e qh - d{q) e - qh ] = (3 - qh)a{q)e qh - (3 + qh)b{q) e - qh (58) 

2{1 + u) [c{q)e qh - d{q)e- qh ] = (1 + 2v - qh)a{q)e qh - (1 + 2v + qh)b{q)e' qh , (59) 

20 



whose solution is 
a(q) = 

b(q) - 
c(q) = 
d(q) - 



2s(q)(l + u) 

2s(q)(l + u) 



\s{q) 



sinh(g/i)e 



-qll 



sinh(2g/i) + 2qh 

smh(qh)e qh 
sinh(2qh) + 2qh 

smh(qh)e~ qh 



\s{q) 



^ + U 'smh(2qh) + 2qh 
smh(qh)e~ qh 



smh(2qh) + 2qh 



(3 + 4z/) 





(60) 




(61) 


sinh(qh)e qh 


(62) 
(63) 


smh(2qh) + 2qh 
sinh(qh)e qh 


smh(2qh) + 2qh 



Solution for a rough bottom 

In the case of a very rough bottom, the relation 
(iv-b), i.e. by: 



2(l + i/) 



c(q)e qh + d{q)e- qh 
The resolution of the four equations then gives 

f-(q) + 2qh 



should be replaced by the condition 
(4z/ - qh)a{q)e qh + (4z/ + qh)b{q) e - qh . (64) 



a(q) 
b(q) 
c(q) 
d(q) 



2s(q)(l + u) 

2s(q)(l + u) 
1- 



f+(q)-2qh 



~\s{q) 



\s{q) 



1 



U(q)f-(q) + M 2 h 2 

(3 + 4v)f-(q) + f + (q)+4qh(l + 2v) 
U(q)f-(q)+^ 2 h 2 

(3 + 4is)f + (q) + f-(q)-4qh(l + 2u) 
f + {q)f-{q)+M 2 W 



where the functions /+ and /_ are defined by 

f ± (q) = l + (3-4u)e ±2qh . 



(65) 
(66) 
(67) 
(68) 

(69) 



Again, as in the 2d case, the standard limit v = ~ can be exceeded without quantitative 
change - except the appearance of negative parts -, up to v = | where the stresses start 
oscillating. 

Semi-infinite medium 

For comparison, in the case of a semi-infinite medium submitted to a ponctual and vertical 
unitary force at r = 0, the stress components are given by Boussinesq and Cerruti's 
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formulae |26|| : 



O, 



<J r 



2vr (r 2 + z 2 fl 2 



1 

2^ 



'1-21/1 



3^r 2 



Z 2 \ 



2^^ l r 2 r 2 (r 2 + z 2 W 2 (r 2 + z 2 ) 3 / 2 

3 rz 2 



a r 



2tt (r 2 + z 2 ) 5 / 2 ' 



(70) 
(71) 
(72) 
(73) 



22 



